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1.  INTRODUCTION 


In  classical  elastoplasticity  (Hill  1950;  Malvern  1969),  it  is  customary  to  assume  the  existence  of  a 
plastic  potential  function  in  stress  space,  which  is  usually  taken  to  be  the  yield  function.  Consistent  with 
this  assumption  is  the  associated  flow  rule,  which  specifies  that  each  component  of  plastic  natural  strain 
increment  is  proportional  to  the  partial  derivative  of  the  potential  function  with  respect  to  the 
corresponding  stress  component.  For  materials  that  work  harden,  the  effect  of  inelastic  straining  on  the 
"size"  of  the  yield  function  must  be  considered.  One  method  of  accomplishing  this  is  to  introduce  the 
concept  of  a-universal  plastic  stress-strain  curve  that  relates  two  scalar  quantities:  an  "effective"  stress  and 
an  accumulation  of  "effective”  plastic  strain  increments  (defined  in  sequel).  This  concept  gives  good 
agreement  with  physical  experiments  for  proportional  loading,  but  less  satisfactory  predictions  for 
combined  stresses  deviating  significantly  from  proportional  loading.  Nevertheless,  this  concept  for 
modeling  work  hardening  is  widely  used,  both  in  analysis  and  in  computational  methodology. 

An  example  of  a  universal  stress-strain  curve  is  the  popular  Johnson-Cook  (1983)  strength  model, 
particularized  to  the  rate-independent,  isothermal  case.*  While  the  parameters  for  this  model  were  derived 
by  fitting  experimental  data,  it  will  be  shown  that  such  a  fit  can  only  be  satisfactory  for  a  limited  range 
of  finite  strains.  An  alternate  approach  is  to  provide  the  universal  stress-strain  curve  as  a  discrete  set  of 
data  points  over  as  wide  a  range  of  effective  plastic  strains  as  desired  and  to  use  an  interpolation  procedure 
for  intermediate  values.  This  report  will  elaborate  on  the  latter  alternative  and  certain  refinements. 

2.  DISCRETE  PAIRS  CONSTITUTIVE  INPUT 

Computer  programs  developed  at  Lawrence  Livermore  National  Laboratory  have  provided  options  to 
input  constitutive  data  as  coordinate  pairs  of  points  on  an  effective  stress  vs.  effective  plastic  strain  curve. 
Specifically,  these  options  apply  to  Material  Types  10  and  24  of  DYNA3D  (Hallquist  1988)  and  Material 
Types  3  and  24  of  NIKE3D  (Hallquist  1984).  An  illustration  of  this  form  of  input  is  shown  in  Figure  1, 
where  experimentally  determined  values  (Weerasooriya  and  Swanson  1991)  for  annealed  OFHC  copper 
are  displayed  as  asterisks.  The  computer  programs  then  calculate  values  of  effective  stress  for  intermediate 
values  of  the  accumulated  effective  plastic  strain  by  linear  interpolation. 


*  The  author  believes  that  any  rate-/temperature-dependant  constitutive  model  should  provide  agreement  with  quasi-static 
experimental  data  when  so  particularized. 
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EFFECTIVE  PLASTIC  STRAIN 

Figure  1.  Effective  stress  vs.  effective  plastic  strain  curves. 
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Figure  1  also  shows  a  corresponding  curve  for  the  Johnson-Cook  model,  based  on  their  published 
parameters  for  OFHC  copper.  The  differences  between  the  two  curves  for  small  to  moderate  effective 
plastic  strains  are  not  regarded  as  significant,  since  there  may  have  been  differences  in  heat  treatments 
between  the  two  sets  of  experiments.  However,  the  Johnson-Cook  model,  which  uses  a  single  power  law 
term,  cannot  represent  the  stress  saturation  that  is  observed  physically.  The  same  limitation  also  applies 
to  the  Zerilli-Annstrong  (1987)  model. 


The  discrete  point  representation  of  constitutive  data  can  also  give  rise  to  undesired  effects.  When 
faced  with  a  limited  number  of  input  points  (16  for  DYNA3D,  8  for  NIKE3D),  one  would  pick  points 
from  available  experimental  data,  not  necessarily  at  uniform  increments  of  effective  plastic  strain,  which 
appear  to  include  the  essential  features  of  these  data.  If  the  experimental  data  are  known  to  represent  a 
smoothly  continuous  phenomenon,  the  first  (and  possibly  the  second)  differences  of  the  selected  input 
points  should  be  checked  for  proper  smoothness. 
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Another  difficulty  can  arise,  at  least  with  the  computational  procedure  for  stress  incrementation 
employed  by  the  DYN A/NIKE  codes.  This  is  illustrated  in  Figure  2,  which  was  obtained  from  an 
NIKE3D  computation  for  the  linearly  increasing  twist  of  a  tube,  using  the  copper  data  from  Figure  1 .  The 
large  "spike"  that  is  observed  at  0.1  s  occurred  when  the  computed  effective  plastic  strain  exceeded  one 
of  the  input  values  (where  there  is  a  slope  discontinuity).  The  problem  is  that  the  constitutive  subroutine 
uses  the  previous  slope  of  the  yield  function  to  evaluate  the  new  yield  stress  and  is  not  aware  of  the  slope 
change  until  the  next  time  step.  In  the  case  of  NIKE3D,  large  time  increments  can  be  employed,  which 
accentuates  die  obvious  discrepancy. 


TIME  (s) 

Figure  2.  "Before  and  After"  stress-time  solution. 

For  the  DYNA3D  code,  where  the  time  increment  is  constrained  to  be  small  for  numerical  stability, 
the  "spike"  also  appears  (perhaps  with  a  smaller  amplitude),  but  is  followed  by  an  oscillatory  response  that 
does  not  rapidly  die  out  This  behavior  has  been  shown  to  be  associated  with  the  occurrence  of  spurious 
negative  plastic  wotk. 
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The  foregoing  phenomenon  can  be  avoided  by  use  of  a  special  algorithm  developed  by  the  author  that 
is  invoked  whenever  the  computed  effective  plastic  strain  reaches  one  of  the  input  data  values. 

3.  THE  "STRADDLE"  ALGORITHM 

The  discussion  of  this  algorithm  refers  to  the  manner  in  which  it  was  implemented  in  the  DYNA3D 
code  for  the  solid  hexahedron  element;  however,  for  other  applications,  the  procedure  is  essentially  the 
same.  The  notation  employed  is  defined  in  the  Nomenclature  section  of  this  report 

Subroutine  F3DM10  takes  the  stress  components  from  the  previous  cycle  and  a  new  set  of  strain 
increments  and  calculates  a  set  of  trial  deviatoric  stress  components  on  the  basis  of  an  elastic  stress 
change.  It  recalls  the  previous  effective  plastic  strain  ep  and  calls  the  inteipolation  subroutine  YIELDS 
to  reconstruct  the  current  yield  stress  Oj  (see  Figure  3).  The  trial  stresses  are  used  to  calculate  the 
von  Mises  effective  (trial)  stress  oT.  If  ctt  <  oj ,  the  stress  change  was  elastic  and  the  trial  deviatoric 
stresses  are  accepted  as  valid.  Otherwise,  the  stress  change  is  plastic  and  the  program  calculates  an 
incremental  effective  plastic  strain  (scalar)  using 


(AeP)* 


T  Y 
3G  +  Hk 


(1) 


In  the  original  version,  of  F3DM10  a  new  yield  stress  is  determined  using 


(o^r  -  oj  +  Hk(Aep)*  (2) 

It  is  this  procedure  that  results  in  a  stress  overshoot  when  ep  and  ep  +  Aep  straddle  a  point  such  as 
ek+1  on  the  yield  curve. 
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EFFECTIVE  PLASTIC  STRAIN 
Figure  3.  Geometry  of  the  straddle  analysis. 


In  the  revised  version  of  F3DM10,  the  new  trial  effective  plastic  strain 

(ej’+1)T  =  e?  +  (Aep)*  (3) 

is  calculated  and  compared  to  ek+1 .  If  >  ek+1,  the  analysis  branches  to  the  straddle  algorithm 

to  deal  with  this  special  case.  It  is  desired  to  determine  a  Aep  and  an  intermediate  slope  He  that  will 

Y 

result  in  a  value  of  Oj+1  which  is  consistent  with  the  input  data.  The  following  conditions  must  be 
satisfied: 


ep 

ej+l 

=  e?  +  Aep 

(4) 

T  _Y 

Aep 

_ a  -  0j 

3G  +  He 

(5) 

Y  Y 

H 

(6) 

e 

ej+l  ej 

Y 

Y  rr  /  P  \ 

(7) 

"  °k+l  +  Hk+l[ej+l  ~  ek+l/ 
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These  simultaneous  equations  may  be  solved  algebraically  to  yield 


H 


e 


and 


Aep  = 


o/-C 


1 


3G  +  Hk+1 


where 


C1  =  °k+l  "  °jY  “  Hk+l(ek+l 


(8) 


(9) 


(10) 


Values  of  ep+1  and  may  then  be  obtained  by  use  of  equations  (4)  and  (7).  Certain  tests  to  verify 
that  £p+1  falls  within  the  range  ek+1  <  e?+1  <  ek+2  are  required. 

Following  these  calculations,  the  normal  flow  of  this  subroutine  is  rejoined  and  the  stress  components 
are  scaled  back  (radial  return  method  [Krieg  and  Key  1976])  so  that  the  new  total  stresses  will  lie  on  the 
enlarged  yield  surface. 

4.  IMPLEMENTATION  OF  THE  ALGORITHM  IN  COMPUTER  CODES 

Although  the  straddle  algorithm  was  formulated  for  use  with  the  solid  element  in  DYNA3D,  it  is  also 
directly  applicable  for  the  NIKE3D  solid  element,  as  well  as  for  the  elements  in  the  2-D  versions  of  these 
codes,  wherever  the  yield  function  is  introduced  as  discrete  points.  The  3-D  programs  also  feature  beam 
and  shell  elements  for  which,  for  certain  options,  discrete  yield  data  are  accepted. 

In  DYNA3D,  both  the  Hughes-Liu  (1981)  and  Belytschko-Tsay  (1981, 1984)  shells  have  this  feature 
and  there  is  a  choice  between  iterative  and  noniterative  plane  stress  plasticity  (Hallquist  and  Whiriey 
1989).  In  NIKE3D,  the  Hughes-Liu  beam  and  shell  foimulations  both  provide  for  discrete  yield  data  and 
both  require  iterative  stress  evaluation  to  satisfy  the  requirement  that  the  stress  component  normal  to  the 
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lamina  system  be  zero.  This  iteration  process  complicates  the  introduction  of  the  straddle  algorithm,  but 
it  has  been  successfully  incorporated  for  the  Hughes-Liu  shell  in  subroutine  s3mns.  In  fact,  the  dashed 
line  solution  in  Figure  2  resulted  from  a  NIKE3D  calculation  after  this  subroutine  was  modified. 

To  date,  the  straddle  algorithm  has  only  been  incorporated  in  those  portions  of  the  DYNA  and  NIKE 
codes  where  a  need  has  arisen.  As  an  example,  the  inclusion  of  this  algorithm  in  DYNA3D  Material 
Model  10  is  listed  in  the  appendix. 

5.  CONCLUDING  REMARKS 

The  straddle  algorithm  was  developed  to  cope  with  nonphysical  "spikes"  and  oscillations  observed  in 
solutions  obtained  using  the  Lagrangian  hydrocodes  NIKE3D  and  DYNA3D,  in  conjunction  with  discrete 
yield  function  input  data.  In  many  applications,  these  perturbations  may  be  only  a  minor  nuisance  in  the 
graphical  output.  However,  when  tracking  material  response  to  transient  loading  in  conjunction  with  a 
material  failure  criterion,  the  decision  as  to  whether  survivability  or  failure  will  result  may  be  altered  by 
the  presence  of  spurious  spikes. 

Consequently,  it  is  argued  that  incorporation  of  the  algorithm  is  an  inexpensive  means  of  avoiding 
unpleasant  surprises  (since  the  number  of  branches  to  the  algorithm  at  each  integration  point  cannot  exceed 
the  number  of  input  data  points  minus  one  regardless  of  the  number  of  cycles  computed). 
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APPENDIX: 

STRADDLE  ALGORITHM  IN  DYNA3D  MATERIAL  MODEL  10 


11 
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The  changes  made  to  DYNA3D  to  incorporate  the  straddle  algorithm  in  the  following  two  subroutines 
are  indicated  in  bold  lower  case  letters. 

SUBROUTINE  F3DM10 (CM) 

COMMON/ BKO 2 /IBURN, DTI, DT2 , ISDO 

COMMON/ AUX2/D1 (128) ,D2 (128) ,D3 (128) ,D4 (128) ,D5 (128) ,D6 (128) , 

1  WZZDT (128) , WYYDT (128 ) , WXXDT (128) ,EINC(128) 

COMMON/AUX11/PO (128) 

COMMON/ AUX1 4/ 

&SIGN1 (128) , SIGN2 (128) , SIGN3 (128) ,SIGN4 (128) , 

&SIGN5 (128) , SIGN 6 (128) , 

&EPX1 (128) , EPX2 (128) , EPX3 (128) , EPX4 (128) ,EPX5(128) ,AUX(128,5) 
COMMON/ AUX1 8 /DD (128 ) , DF (128 ) 

COMMON/ AUX1 9 /SP , BFAC (128) ,DR1V(128) ,DR2V(128) ,W1 
COMMON/ AUX2  0 / 

&  AJ2 (128) ,SJ2 (128) ,SCALE(128)  , FJL (128)  , CC (128 ) ,P (128) , 

&  POLD (128) ,DAVG (128) ,FJK(128) ,SPECEN(128) ,AK(128) ,AKT(128) , 

&  YWH (128) 

COMMON/AUX33/IX1 (128) , 1X2 (128) , 1X3 (128) , 1X4 (128) , 1X5 (128) , 

1  1X6(128) ,1X7 (128) ,1X8 (128) ,MXT(128) , NMEL 

COMMON/ AUX35/RHOA (128 ), CXXA (128 ) 

COMMON / AUX3  6 / LFT , LLT 
COMMON/EOSD/PC (128) ,SHRM(128) 
logical  plastic (128) 
common/ state/plastic 
DIMENSION  CM(1) ,il (128) 

DATA  THIRD/. 3333333333333/ 

MX=  48* (MXT (LFT) -1) 

SP=0 . 0 
G=CM (MX+1 ) 

QH=CM (MX+3) 

QS=CM (MX+2 ) 

A1=CM (MX+5) 

A2=CM(MX+6) 

ISPALL=CM (MX+7 ) -1 . 

NPT=CM (MX+1 6 ) 

G2=2 . *DT1*G 
DO  10  I=LFT, LLT 
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10  DAVG (I) “-THIRD* (D1 (I) +D2 (I) +D3 (I) ) 

DO  20  I=LFT, LLT 
SHRM(I) =G 
AKT ( I ) =QH 
PC (I) =EPX5 (I) 

20  AK(I) =QS+QH*EPX1 (I)  +  (Al+A2*PO (I) ) *MAX  (0 . , PO  (I) ) 

DO  30  I=LFT, LLT 

SIGN1 (I) =SIGNl (I) +PO (I) +G2* (D1 (I) +DAVG (I) ) 

SIGN2 (I) =SIGN2 (I) +PO (I) +G2* (D2 (I) +DAVG (I) ) 

SIGN3 (I)=SIGN3 (I) +PO(I) +G2* (D3 (I) +DAVG (I) ) 

SIGN 4 (I) =SIGN4 (I)  +  . 50*G2*D4 (I) 

SIGN5 (I) =SIGN5 (I) +. 50*G2*D5 (I) 

30  SIGN 6 (I) =SIGN6 (I) +. 50*G2*D6 (I) 

IF  (ISPALL.NE.2)  GO  TO  80 
DO  40  I=LFT, LLT 
40  DAVG ( I ) =PO ( I ) -EPX5 ( I ) 

DO  50  I=LFT , LLT 

50  SCALE (I) = . 50* (1 . +SIGN (1 . , DAVG (I) ) ) 

DO  60  I=LFT, LLT 

PC ( I ) “SCALE ( I ) *EPX5 ( I ) 

SIGN1 (I) =SCALE (I) *SIGNl (I) 

SIGN2 (I) =SCALE (I) *SIGN2 (I) 

SIGN3 (I) “SCALE (I) *SIGN3 (I) 

SIGN4 (I) “SCALE (I) *SIGN4 (I) 

SIGN5 ( I ) “SCALE ( I ) *SIGN5 (I) 

60  SIGN6 (I)=SCALE (I) *SIGN6 (I) 

DO  70  I=LFT, LLT 
70  EPX5 (I) =PC (I) 

80  IF  (NPT.EQ.0)  GO  TO  100 

CALL  YIELDS (NPT, CM (MX+17) , CM(MX+33) ,EPX1, AK, AKT, LFT,  LLT,  il) 

100  DO  110  I=LFT, LLT 

AJ2 (I)=. 5* (SIGN1 (I) **2+SIGN2 (I) **2+SIGN3 (I) **2) +SIGN4 (I) **2+SIGN5 
1  (I) **2+SIGN6 (I) **2 
S J2 (I) =SQRT (3 . *AJ2 ( I) ) 

110  CONTINUE 

do  115  i  =  1ft, lit 
plastic (i) =  . false . 

if (ak(i) .le.sj2 (i))  plastic (i)  =  .true. 

115  continue 
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C....  CALCULATE  INCREMENT  IN  EFFECTIVE  PLASTIC  STRAIN  &  NEW  YIELD  STRESS 
THRG=THIRD/G 
DO  120  I=LFT, LLT 

if (plastic (i) )  then 

CC (I) =DIM (SJ2 (I) , AK (I) ) / (1 . E-30+DIM(l . , -THRG*AKT (I) ) ) 

CC (I) =THRG*AMIN1 (SJ2 (I) ,CC(I) ) 

epxt=epxl (i) +cc (i) 
anc-fil  (i)+16 

if (epxt .gt.cm(32) )  go  to  120 

if  (epxt.gt.cm(m)  .and.  epxt . le . cm (m+1 ) )  then 
hb=  (cm (m+1 7)  —cm  (m+1 6)  )  /  (cm (m+1)  -cm (m) ) 
cl=cm (m+1 6)  -ak  (i)  -hb*  (cm(m)  -epxl  (i) ) 
akt  (i)  =  (3 .  *g*cl+hb*  (s  j2  (i)  -ak  (i) ) )  /  (s  j2  (i)  -ak  (i)  -cl) 
cc  (i)  =  (s  j2  (i)  -ak  (i)  -cl)  /  (3 .  *g+hb) 
ak  (i)^nax  (cm(m+16)+hb*  (epxl  (i)+cc(i)  -cm(m) ) ,  0 . ) 
elseif (epxt . gt . cm (m+1 ) )  then 
write (13, 1)  i 
call  adios(2) 
else 

AK ( I ) =MAX ( AK ( I ) +AKT ( I ) *CC  (I ) , 0 . ) 

endif 

else 

cc(i)s  0. 

endif 

120  CONTINUE 

C....  CALCULATE  FACTOR  TO  SCALE  STRESSES  BACK  TO  YIELD  SURFACE 
SEP  S= 1 . E-3  0  *G 
DO  130  I=LFT, LLT 

SCALE (I) = (AK (I) +SEPS) / (MAX (AK (I) , SJ2 (I) ) +SEPS) 

130  CONTINUE 

C....  SCALE  BACK  STRESS  AND  INCREMENT  PLASTIC  STRAIN 
DO  140  I=LFT, LLT 
SIGN1 (I) =SCALE (I) *SIGN1 (I) 

SIGN2 (I) =SCALE (I) *SIGN2 (I) 

SIGN3 (I ) =SCALE ( I ) *S IGN3 ( I ) 

SIGN4 (I) -SCALE (I) *SIGN4 (I) 

SIGN5 ( I ) =SCALE (I) *S IGN5 ( I ) 

SIGN 6 (I) —SCALE (I) *SIGN6 (I) 
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EPX1 (I) =EPX1 (I) +CC (I) 

140  CONTINUE 

IF (ISPALL.NE. 1)  RETURN 
DO  150  I=LFT,LLT 

AJ2 (I) =. 5* (SIGN1 (I) **2+SIGN2 (I) **2+SIGN3 (I) **2) +SIGN4 (I) **2+SIGN5 
1  (I)**2+SIGN6(I)**2+1.E— 12 

150  SJ2 (I) =SIGN1 (I) *SIGN5  (I) **2+SIGN2 (I) *SIGN6 (I) **2+SIGN3 (I) *SIGN4 (I) 
1  **2-SIGNl (I) *SIGN2 (I) *SIGN3 (I) -2. *SIGN4 (I) *SIGN5 (I) *SIGN6 (I) 

DO  160  I=LFT, LLT 

AKT (I) =-SQRT (27 . /AJ2 (I) ) *SJ2 (I) *0 . 5/AJ2 (I) 

AKT (I)=SIGN (AMIN1 (ABS (AKT (I) ) , 1 . ) ,  AKT (I) ) 

160  YWH ( I ) =ACOS (AKT ( I ) ) *THIRD 
DO  170  I=LFT, LLT 

170  SJ2 (I)=2.*SQRT(AJ2 (I) *THIRD) *COS (YWH (I)  ) 

DO  180  I=LFT, LLT 
180  DAVG ( I) =PO (I) -S J2 (I) -EPX5 (I) 

DO  190  I=LFT, LLT 

190  SCALE (I ) = . 50* (1 .+SIGN (1 . , DAVG (I) ) ) 

DO  200  I=LFT, LLT 

PC (I)  =SCALE (I) *EPX5 (I) 

SIGN1 (I) =SCALE (I) *SIGN1 (I) 

SIGN2 (I) =SCALE (I) *SIGN2 (I) 

SIGN3 (I) =SCALE (I) *SIGN3 (I) 

SIGN4 (I) =SCALE (I) *SIGN4 (I) 

SIGN5 (I) =SCALE (I) *SIGN5 (I) 

200  SIGN6 (I) “SCALE (I) *SIGN6 (I) 

DO  210  I=LFT, LLT 
210  EPX5 (I) =PC (I) 

RETURN 

1  format ('  EXCESSIVE  PLASTIC  STRAIN  IN  ELEMENT' , ilO , / ) 

END 

SUBROUTINE  YIELDS (NPT, EPS, SIGE, EPX, AK, QH,  LFT,  LLT, il) 

C  Rewritten  by  Glenn  Randers-Pehrson,  BRL,  6  April  1992,  to 

C  interpolate  correctly  within  the  beginning  and  end  segments, 

C  to  extrapolate  correctly,  and  to  run  faster. 

C 
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c 


INTEGER 

REAL 


NPT, LFT, LLT 

EPS (16) , SIGE (16) , EPX (LLT) , AK (LLT) ,  QH (LLT) 


INTEGER  I,  L, IL (128) 

REAL  SLOPE (16) 

DO  10  I=LFT, LLT 
IL  (I) =0 

IF (EPX (I) .LT.EPS (2) )  IL(I)=2 
IF (EPX(I) . GT . EPS (NPT— 1 ) )  IL(I)=NPT 
10  CONTINUE 

DO  30  L=2, NPT 
DO  20  I=LFT , LLT 

IF  (IL (I)  .EQ.O .AND.  (EPX (I)  .GE.EPS(L-l)  .AND. EPX (I)  .LE.EPS (L) ) ) 
+  IL (I)=L 

20  CONTINUE 

SLOPE (L)  =  (SIGE (L) -SIGE (L-l)  )  /  (EPS (L) -EPS (L-l) ) 

30  CONTINUE 

DO  40  I=LFT, LLT 
QH (I) =SLOPE (IL (I) ) 

AK (I) =SIGE (IL (I) — 1) +QH (I)  *  (EPX (I) -EPS (IL (I) -1) ) 

40  CONTINUE 

RETURN 

END 
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NOMENCLATURE 


H„ 


'm 


()T 

()* 


-  °k+l  ~  °J  ~  Hk*l(ek 


.»  -  «r) 


shear  modulus 

slope  determined  by  the  straddle  algorithm 


Hk  =  (°k+i  “  °k)/(£k+t  "  ek);  of  yield  curve  for  interval  ^  to  ekH 

Hk+1  =  (°k+2  "  °k+l)/(^k+2  “  ^k+l) 

N  number  of  input  (em,  om)  pairs 

Aep  incremental  change  in  the  effective  plastic  strain 

input  values  of  effective  plastic  strain,  m  =  1,  2, ....  N 


—AefjAe?: ;  current  value  of  effective  plastic  strain 
N|  3  J  •> 


=  E 

current  value  of  yield  stress 


input  values  of  effective  stress,  m  =  1,2,  •••,  N 


3  T  T 

,  _Tjj  Tjj  ;  trial  value  of  von  Mises  effective  stress 


deviatoric  components  of  Cauchy  stress  tensor 
trial  value 

quantity  calculated  using  original  program 
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